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We investigate the 2D instability recently discussed by Gallet et al. (2010) and Ilin 
& Morgulis (2013) which arises when a radial crossflow is imposed on a centrifugally- 
stable swirling flow. By finding a simpler rectilinear example of the instability - a sheared 
half plane, the minimal ingredients for the instability are identified and the destabiliz¬ 
ing/stabilizing effect of inflow/outflow boundaries clarified. The instability - christened 
‘boundary inflow instability’ here - is of critical layer ty pe where this layer is ei ther at 
the inflow wall and the growth rate is O(y^) (as found bv lllin fc Morgulis I (|2013l) ). or in 
the interior of the flow and the growth rate is 0{r] log l/rj) where r] measures the (small) 
inflow-to-tangential-flow ratio. The instability is robust to changes in the rotation profile 
even to those which are very Rayleigh-stable and the addition of further physics such as 
viscosity, 3-dimensionality and compressibility but is sensitive to the boundary condition 
imposed on the tangential velocity field at the inflow boundary. Providing the vorticity is 
not fixed at the inflow boundary, the instability seems generic and operates by the inflow 
advecting vorticity present at the boundary across the interior shear. Both the primary 
bifurcation to 2D states and secondary bifurcations to 3D states are found to be super¬ 
critical. Assuming an accretion flow driven by molecular viscosity only so ry = 0{Re~^), 
the instability is not immediately relevant for accretion disks since the critical thresh¬ 
old is and the inflow boundary conditions are more likely to be stress-free 

than non-slip. However, the analysis presented here does highlight the potential for mass 
entering a disk to disrupt the orbiting flow if this mass flux possesses vorticity. 


1. Introduction 

Rotating flows are ubiquitous in nature and industrial applications so understanding 
their stability continues to be an important and active area of research. The extent of this 
activi ty is perhaps best epitomised by the huge bo dy of work studying Taylor-Couette 
flow ( Couettel llBSSt Mallock 1888t Tavlor 1923ll - the flow between two concentric 
cylinders rot ating at different r a tes - which has beco me the laboratory paradigm of the 
subject (e.g. Tuckerman I ( 20141) : IPardin et~^ ( 2014 1 and references herein). Despite all 
this work, however, only very recently has it been realised that imposed radial flow ca n 
destabilise an otherwise stable r otatin g flow ( Gallet et al ] l201(ll : lllin fc Morerulis iboioll . 
The paper bv lllin fc Morguhsi ( 2013h is particularly revealing because it demonstrates 
that the (non-dimensionalised) flow 

U(r,6») = 

r r 

































2 


R. R. Kerswell 


between two concentric cylinders at r = 1 and r = a > 1 is inviscidly unstable to 
infinitesimal 2D oscillatory disturbances for either inflow 77 > 0 or outflow 77 < 0 providing 
77 is not too large. This is surprising because, firstly, in the absence of radial flow, this 
rotating flow is (marginally) ‘Rayleigh-stable’ as the angular momentum of the flow 
nowhere decr eases in magnitude radially (note this is s trictly a condition on axisymmetric 
disturbances: Ravleighi (ll917ll : lDrazin fc Reid ( 198lh §15.2). Secondly the flow also fails 
a require ment for 2D inst a bility - the rotating flo w version of Rayleigh’s inflexion point 
theorem ( Rayleigh ( 188(t) : Drazin &: Reid ( 198lll §15.3 and problem 3.2 p 121). Thirdly, 
it is also somewhat counterintuitive that the instability occurs for both converging and 
diverging radial flows since the former are stable an d the latter ge nerally unstable when 
the rotation is a bsent as in Jeffery-Ham el flow (e.g. IPrazin I (Il999l) b 

The results of[li in fc Morgiilis I (|2013h also possess many intriguing features of which 
four stand out. Firstly, the existence of an imposed normal flow through the cylinder 
walls increases the order of the linear operator describing the evolution of small inviscid 
disturbances from the normal 2 to 3. This means that an extra boundary condition has 
to be imposed beyond the usual no-normal-flow conditions at either cylinder wall. While 
it is straig htforward to argue th at this extra condition must be imposed at the inflow 
boundary ( Ilin fc Morgulis ^2013^ . predicting the effect of a specific choice is less clear: if 
a no-slip condition is imposed there is instability whereas a vanishing v orticity condition 
gives stability. Secondly, in the limit of vanishing radial flow (77 —>■ 01. iRin &: Morguli^ 
( 2 OI 3 II find growth rates which scale as 0{^) rather than the generic 0(77) one might 
expect. This unusual growth rate scaling arises because there are no discrete modes of 
the linear problem which can satisfy the usual 2 no-normal velocity conditions for 77 = 0 . 
This gives rise to a non-standard singular perturbation analysis, the robustness of which 
is unclear to, say, changes in the rotation profile and/or to the addition of extra physics. 
Thirdly, the 77 —» 0 asy mptotics presented seems incomplete. The instability described in 
iRin fc Morgulis I ( 2013l l has a critical layer character but only where this critical layer is 
actually at the inflow boundary. This suggests the existence of further instabilities with 
an interior critical layer separated from the boundary. Lastly, the mechanism for the 
instability is unclear. Crossflow and shear would seem obvious ingredients but rotation 
or curvature not necessarily so. It is also not apparent whether the energy to feed the 
instability comes wholly ’through the boundary’ or is extracted at least in part from the 
(interior) rotational energy of the underlying flow. 

To keep this study manageable, the focus here is on the 0 < 77 <g; 1 situation which 
repre sents small radial in flow on a predominantly rotating flow. The motivation for this 
(as in Gallet et al.l ( 201^ 1 is the accretion disk problem where certainly for cold and hence 
wea kly-ionised disks, t he source of inferred d i sk turbul e nce remains a hotly contest ed issue 
(e.g. Dubrulle et al. ( 2005 ): Shariff ( 2009| ): iBalbus I ( 201 lh : [ji fc Balbiisl ( 201,3h ). As a 
result there is considerable interest in uncovering robust linear instability mechanisms. 
Interestingly, the existence of the radial accretion flow is rarely included in theoretical 
models since it is so small - 0(l/i?e) sma ller than the azi muthal flow where Re is huge 
when based on a molecular viscosity (e.g. IPubrullel ( 1992 ) quote figures of 10 ^"^ — 10 ^®) 
- and presumab ly its presence only felt over 0{Re) timescales which are far too large to 
be relevant (e.g Shariff ( 2009h estimates that it would take longer than the age of the 
universe f or molecular viscosity to diffuse momentum across a typical disk). However, the 
results of lllin &: Morguli^ ( 20131) suggests that such a flow could actually drive linear 
instabilities over a much shorter 0{\/Re) timescale. 

The plan of the paper is to start with perhaps the simplest example of the instability 
which is just a sheared half plane of fluid with imposed inflow. The effect of adding 
viscosity is discussed as well as the introduction of an outflow (suction) boundary so 
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that the flow domain becomes a channel. Then the discussion turns to rotating flow with 
more general profiles, including the solidly Rayleigh-stable Keplerian profile U = Ify/rO, 
to examine the robustness of the instability. The effect of further physics in the form of 
viscosity, 3-dimensionality and compressibility are also broached. Finally, the nonlinear 
aspects of the instability are probed ranging from a weakly nonlinear analysis around 
the primary 2D bifurcation through to secondary bifurcations and the ensuing fully 3D 
finite amplitude solutions. 

The findings of the paper, organised under the various questions posed above, are as 
follows. 

(a) Is curvature or rotation important for the instability?. No, all the features of the in¬ 
stability are reproduced in a rectilinear flow with inflow described in ^2.11 The instability 
operates by advecting a source of vorticity at the inflow boundary across shear. 

(b) Are there other (non-boundary-layer) modes of instability caused by an inflow 
boundary? Yes, instabilities exist with interior critical layers distinct from the inflow 
boundary: see ^2.1.21 ^3.21 and expressions (12.121) . (I2.48|) and (13.261) . Their growth rates 
are 0{r] \og{l/r]) ), w hile much larger than O in). are smaller than the inviscid boundary- 
layer modes fou nd by Ilin fc: Morgulis ( 2013ll and the equivalent viscous boundary layer 
modes found bv iGallet et al. ( 2010 ). 

(c) Given the sensitivity to what is chosen for the extra boundary condition, how 
generic is the instability across possible boundary conditions ? For finite 77 the instability 
looks generic with only a no-vorticity boundary condition obviously ensuring stability: 
see mi However for 0 < ry <C 1, any restriction on the normal derivative of the tan¬ 
gential velocity effectively kills the instability: see )12.1.4I The non-slip condi tion always 
seems to allow instability to occur ( Gallet et al 1 120101: lllin fc Morgulis 112011 . 

(d) How robust is the instability to different rotation profiles? Yevy robust. The form of 
the shear is unimportant for the boundary-layer instability and of secondary importance 
for the critical-layer instability: see T3.21 

(e) What is the effect of adding viscosity? The presence of viscosity introduces a thresh¬ 
old crossflow of 0(Re~^) in the rectil i near s ituation where long streamwise wavelengths 


are permissible ( Nicoud fc AngilelU ( 1997 ) and ( 12 . 21 ) o r 0(Re f or the rotational 


situation where the azimuthal wavenumber is an integer ( Gallet et al.l ( 2010 l) and 

(/) What is the effect of adding 3 dimensionality? Squire’s Theorem effectively holds 
for the boundary layer instabilities since only the shear at the boundary is important 
and curvature is secondary. As a result, adding 3 dimensionality leads to less unstable 
disturbances: see in 

{g) How robust is the instability to adding further physics? The instability survives the 
addition of compressibility (see S]) , 3-dimensionality (see (l5|) and viscosity (see (12.21 and 
|). It is also insensitive to the exact shear present as long as it is non-vanishing (see 


(/i) Has this instability bee n seen before Gallet e t all \201(\ ) and, Ilin & Moraulis I 
Yes, in a rectilinear form bv iNicoud fc Angilella ( 19971) who studied ‘generalised plane 
Couette flow’ where a streamwise pressure gradient is imposed to counterbalance the 
effects of crossflow in the streamwise momentum balance (see (J2]). Doering et al. I ( 200r][) 


also saw the instability without an imposed pressure gradient. In this case the intro¬ 
duction of crossflow not only adds a new flow component to the base state but also 
changes its cross-stream shear. From a stability perspective, how these two effects in¬ 
teract can be subtle and ‘suction’ (as it is t ypically called ) can be found to stabilise 
or des t abilise existing shear ins t abilities (e.g. Hai ns I ( 19711) : iFransson and Alfredsson '' 


( 2003h : Guha fc Frigaard] ( 2011)1 ): Deguchi et al.l ( 2014 ) ). The situation is similar in the 


Taylor-Gouette problem where radial flow can either stabilise or destabilise the well 
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known Tavlor vortex instabilitv (IChang & Sartorv 1967t Bahl 1970 

: Min & Lnentow 

1994 Kolvshkin Kr, Vaillancourt Il997l: iKolesov Shauakidze 1 199 


Serre et al. 

I 2 OO 8 : 

Martinand et al. 120091) . The imoortance of the work of Gallet et al. ( 

2010t) and Ilin & Morgulis 

( 2013|) is that thev studied the effect of radial flow on centrifugallv-stable Tavlor-Couette 


flow. 


(i) Is the instability supercritical or subcritical? It is a supercritical instability for 
bifurcations where the crossflow is fixed and non-slip boundary conditions are applied at 
the boundaries in the presence of viscosity: see 116.11 and Appendix B. 

(j) Can secondary bifurcations off the 2D solutions reach crossflow values which are 
below the threshold for (primary) instability? There is no evidence for this. The six 3D 
bifurcations found are all supercritical leading to even higher crossflow values. 

(fc) What is the energy source for the inviscid instability? The instability draws its 
energy from the underlying shear. The energy of this is replenished by the pressure field, 
which drives the crossflow, doing work. 

(Z) Is this instability possibly relevant for accretion disks? In a quiescent disk, the 
(molecular) viscosity-driven accretion flow is 0{Re~^) whereas the critical threshold for 
linear instability is a radial flow of 0{Re~^/^) leaving aside any issues about the exact 
form of the outermost boundary conditions. This, together with the fact that no signs 
of subcriticality have been found in either the primary instability or of any secondary 
bifurcations, indicates that the instability is not operative in isolation . However, if some 
other process is able to generate a larger radial flow, then this instability may be trig¬ 
gered as a secondary consequence. The instability primarily derives its energy from the 
gravitational body force working on the accreting flow. 


Ilin and Morgulis have themselv es continued their inve stigation to consider th e effect 
of viscosity ( Ilin fc Morgulis 20151) and 3-dimensionality ( Ilin &: Morgulis 20I5tJ) albeit 
from a complementary perspective: given a crossflow what is the smallest critical swirling 
flow for instability? This is equivalent upon rescaling to the problem of fixing the swirling 
flow and finding when instability disappears as the crossflow is increased rather than 
decreased as studied here. That there is a finite range of crossflow for instability when 
there is both an inflow and outflow boundary appears a generic observation. There is 
also a tempting general interpretation - the inflow boundary is responsible for initially 
destabilising the flow as the crossflow is increased in magnitude from zero but ultimately 
the outflow (or more commonly labelled the ’suction’) boundary stabilises the flow again 
when the crossflow becomes large enough. The particularly simple half-plane problem 
treated below helps m otivate this simple characte r isation. We hencefort h refer to the 
instability first seen by Nicoud &: Angilella I ( 1997ll . [Poering et al. I ( 2000ll . Gallet et al 
( 20inll and lllin fc Morgulis I (l20l¥ as the ‘boundary inflow instability’. 


2. 2D Instability in Simplest Form: The Half Plane 

The boundary inflow instability operates by advecting vorticity acr oss shear, occurs 
in 2D and does not need viscosity or underlying vorticity as shown in Ilin &: Morgulis 1 
( 2m3h . Here, we demonstrate that curvature or rotation are not necessary either by 
discussing a rectilinear example of the instability. The instability does need a boundary, 
however, so cannot be captured by a local analysis. To see this, consider the simplest 
possible set-up: a 2D shear flow 


U = yx + py 


( 2 . 1 ) 
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where there is a constant pressure gradient —77X and U conveniently solves both Euler 
equations and the Navier-Stokes equations (by desigiij). The linearised inviscid equation 
for the perturbation vorticity lo = —A'i/;, where if is the streamfunction (u = V x ifz), is 


/ d d d\ ^ 

which represents just advection of the vorticity and has solution 

w(a;, y, t) = ujo{x - yt + - yt) 


( 2 . 2 ) 


(2.3) 


where ojQ{x^y) is the initial vorticity distribution. In an unbounded domain, wq can be 
Fourier transformed so that it is sufficient to consider only ujq = exp(i(fcx + ly)). Then 
the stream function is 


t-2 I 72 1 „ 

If 4.\ /i “T t i(kx-\-(l — kt)y-\-r}\T:kt'^—It]) 


(2.4) 


which exhibits tran sient growth but no linear instability just as in the rj = 0 case (jOrr 


19071: iFarrell 11198711 . 


2.1. Half Plane: Inviscid 

A half plane, however, can permit a starting vorticity distribution which spatially grows 
towards the boundary. Consider a boundary at y = 0 and 77 > 0 so that this is an inflow 
boundary when the fluid domain is y > 0. Taking the Fourier transform in x of (lOl) 
forces 

A7/>(x, y, t) = f{y - (2.5) 

where / is an arbitrary function. If ip is to be a modal disturbance which depends 
exponentially on t (and x), i.e. ifix, y, t) = ipkiv) exp{ikx + at), then 

/(^) = ( 2 . 6 ) 


is the only possibility {A is an arbitrary normalisation) and therefore 


dy'^ 


(2.7) 


(no modal solution exists for 77 = 0). We will show that this equation, where the initial 
vortical distribution decays exponentially as y —> 00 , has growing modal disturbances 
(i.e. 5fte((T) > 0). No instability is possible if the boundary is an outflow boundary - i.e. 
?7 < 0 - as A is forced to be 0 by boundedness. This suggests that the instability exists 
providing the boundary is not a zero-vorticity boun dary which would force A = 0. In what 
follows, we adopt a non-slip boundary condition as lllin fc Morgulis I (l2013^ originally did 
but will revisit this issue in T2.1.4I (equation (12.71) is a third order differential equation 
for Tp integrated once to include an arbitrary constant - hence 3 boundary conditions are 
needed to specify a unique solution). 

To confirm there is instability, must be solved to derive the dispersion relation. 
Using a Green’s function approach, setting A to 1 (w.l.o.g.) and imposing the 2 (usual 
inviscid) boundary conditions that = 0 and limy_>oo V’fc(?/) = Oj (|2.7|) has the 


f Removing the pressure gradient means that the base state is no longer a constant shear in 
the viscous situation. 

















6 


R. R. Kerswell 


solution 


Jo ^ Jy ^ 

( 2 . 8 ) 

A third boundary condition needs to be applied to give the dispersion relation, and with 
no-slip at the influx boundary {d-ipk/dy = 0 at j/ = 0 ), this is the relation 


f 


a-^y- 


{av+^iky^)/r^^- ^ 


(2.9) 


There is no solution to this for fc = 0 indicating the absence of a ID instability but 
there is instability in 2D (fc 7 ^ 0). The expression (12.91) is valid for finite y but it is now 
useful to consider small y where this instability can be understood as a critical layer 
instability. The (special) case where this critic al layer is at the ( influx ) boundary (so it 
is in fact a boundary layer) was discussed by Ilin fc Morgulisi ( 20131) in their inviscid 
Taylor-Couette set-up. The growth rates for such modes are 5ie(cr) = 0{yjy). The other 
(generic) case when the critical layer is distinct from the inflow boundary has not been 
discussed before. The growth rate here is smaller - 5ie(cr) = 0(77 log(l/? 7 )) (see below) - 
but is still larger than the default 0 ( 77 ) which would be expected. 

The asymptotic form of the dispersion relation (12.91) can be derived as 77 —> 0'^ using 
standard steepest descent/saddle point ideas. Here we need two parts of the integrand to 
contribute at the same leading order and precisely cancel. This can happen in two ways 
since there is just one saddle point at ys = ia/k: the contribution from an end point 
(clearly y = 0) cancels the contribution from the saddle point (the interior critical layer 
case) or the saddle point and the end point are effectively one and the same asymptotically 
(the boundary layer case). In the former case, the leading contributions from the end point 
at y = 0 (1st term on LHS in (12.10(1 ) and from the saddle point (2nd term on LHS) must 
satisfy 


ik 


( 2 . 10 ) 


Now, for the saddle point to be asymptotically separated from the endpoint y = 0, |cr/A:| = 
0(1). This together with the fact that the magnitude of the saddle point contribution 
0{y^ exp{ar(Ji/{kr])) (where a = Ur + i(Ji) has to be 0 ( 77 ) requires either Or = 0(1) and 
then (Ti <C 1 or (Tr 1 and at = 0(1). The former case is inconsistent because the saddle 
point is in the wrong part of the y—complex plane leaving the contribution from the end 
point y = 0 unbalanced. The latter situation, however, does yield solutions. Defining 

a = iai + S{ri)dr ( 2 - 11 ) 

with (Ti and dr both 0 ( 1 ) quantities and ^( 77 ) —>■ 0 as 77 —>■ 0 , then (| 2 . 10 l) requires 


5{ri)dr =flog-+ log ) -^ 77 -^ 0 ( 77 ^ log 1 / 77 ) 

2ai\ r] k ' 


and 


Arg 


J{a^/{2kr])+Tv/4) 


= 0(77 log 1 / 77 ) 


( 2 . 12 ) 


(2.13) 


so there is instability with asymptotic growth rate 0(77 log I/ 77 ) for a discrete set of 
frequencies 


= —'\J ^Trkr]{8n — 1 ) 


(2.14) 
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asymptotic 


actual 


n 

hi 

(Tt* 

hi 

(T f 

100 

-50.1012 

0.17909 

-50.10143 

0.1790934 

10 

-15.7539 

0.42268 

-15.75720 

0.4226325 

9 

-14.9387 

0.43871 

-14.93867 

0.4386500 

8 

-14.0684 

0.45724 

-14.07266 

0.4571595 

7 

-13.1449 

0.47903 

-13.14981 

0.4789312 

6 

-12.1513 

0.50526 

-12.15722 

0.5051295 

5 

-11.0690 

0.53781 

-11.07620 

0.5376254 

4 

-9.8686 

0.57997 

-9.877886 

0.5796846 

3 

-8.5004 

0.63820 

-8.513168 

0.6377301 

2 

-6.8647 

0.72800 

-6.884666 

0.7270530 

1 

-4.6894 

0.90317 

-4.732350 

0.9003686 


Table 1. The 10 most unstable eigenvalues from the dispersion relation (12.151) and asymptotic 

estimates from (12.1711 . 


where n is an integer of 0 ( 1 / 77 ). The form of the corresponding eigenfunction is discussed 
in ^TT 21 

The other situation - the boundary layer case - arises when the saddle point is within 
0{^Jrf) of the endpoint at 7 / = 0, i.e. \a\ = O(y^). At this point, the endpoint and the 
saddle point contributions cannot be considered separately. Instead y and a must be 
rescaled so that (1^ becomes 

e-^^-^^^dz = 0 (2.15) 

where z := y/\j2ri/k and 

(J = yjkri/2 (dr + ia-i) + 0 ( 77 ). (2-16) 

This must be solved numerically but a good estimate for the eigenvalues can be found by 
treating the contributions from the end point and saddle point as if they are separated. 
This means taking the integer tt, to be 1 <C n <C I/77 in the frequency expression (12.1411 for 
the internal critical layer mode and calculating the corresponding growth using ( 12 . 1211 . 
This leads to the asymptotic form 

(Ti ^ —y/'jT{8n — 1) ar ~ log(7r(T^) as integer —>• 00 (2-17) 

Wi\ 

which performs very well even for the first eigenvalue since di is already < —4.7 then: 
see Table 1. Figure [T] shows a typical critical layer eigenfunction and 3 boundary layer 
eigenfunctions for fc = 1 and 77 = 10 “^. 



2.1.1. Inviscid Asymptotics 0 < 77 <g; 1 for the Boundary Layer Instability 


Here a modal streamfunction solution of the form '^tjx^ y.t) = 'i/)( 7 /)e^^^~*~° ’ * is so ught 
with a boundary layer of thickness y/rj as first uncovered by lllin fc Morgulis I (j2013t l (see 


their §3.2). Defining the new boundary layer variable f := y^Jk/{2r]) and splitting the 
streamfunction into an expansion of large-scale parts and a boundary layer corrections 
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Figure 1. Unstable boundary layer eigenfunctions (upper plot: blue solid line is the most 
unstable, the red dashed line is the third most unstable and the black dash-dot line is the fifth 
most unstable mode) and an unstable critical layer eigenfunction (lower plot with the critical 
layer shown as a red dashed line) for 77 = 10“® and fc = 1. In all, the real part oi u = ipy is 
shown. 


(hatted variables) 

V' := {'4’oiy) + i’oiO ) + Vvifpiiv) + 'fpiiO ) + ■•■, (2-18) 

we look for an instability with vanishing frequency (=fcx speed of the inflow boundary) 
to leading order, 

a := ^/r]kj2a + 0{ri). (2-19) 

The governing equation 

{a + iky + ridy){'tjjyy - k'^'ijj) = 0 ( 2 . 20 ) 


is then simplified to 


i’Oyy - k'^4’0 = 0 


(2.21) 
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for the leading flow and 

(a + 2i^ + = 0 (2.22) 

for its boundary layer correction. The interesting observation here is there is no large- 
scale solution -fo which can handle both the rj = 0 boundary conditions that v{x, 0, t) = 0 
and \irny^oov{x,y,t) = 0. Crucially, this means that the boundary layer correction ifo 
must be 0(1) so as to contribute at leading order to fix up the v{x,0,t) = 0 boundary 
condition rather than the usual 0{y/rj) to satisfy the extra rj 0 tangential boundary 
condition u{x,0,t) = 0. As a consequence, there is an 0(1/^) tangential velocity u 
in the boundary layer which must vanish at ^ = 0. Since the large-scale flow cannot 
contribute at this order, this non-slip condition is solely on the boundary layer flow and 
is sufficient to determine the growth rate. Integrating ()2.22|1 twice and incorporating the 
fact that lim^_>oo ' 00 ? = 0, leads to 

pOO 

= (2.23) 

where A is an arbitrary constant. Imposing it = 0 at y = 0 then forces '0o5(O) = 0 which 
is precisely condition (I2.15|) . 


2.1.2. Inviscid Asymptotics 0 < ry -C 1 for the Critical Layer Instability 

An interior critical layer instability is constructed by looking for a mode of 0(1) fre¬ 
quency. Adopting the expansion (I2.11|) where again both ai and ar are 0(1), the critical 
layer is centred on y^ ■= —di/k (so Ui < 0) and as in the boundary layer case, has 
thickness 0{^). Defining the critical layer variable 


. ^ y-Vc 

Vv 


(2.24) 


the (leading order) streamfunction in the critical layer, 0, satisfies 


{Sar/^/q + ikf, + = 0 (2.25) 

which can be integrated once to give 

= ^e-^^A/An-^ke/2 ( 2 . 26 ) 

(choosing the normalisation of the mode here for convenience later) and then two further 
times to give 


4 ,^ = dx + A^ ^^(^±^l^ + A^+h.o.t. asC 


y ±00 
(2.27) 


if = y/q 



+ B 





A 


+ B + h.o.t. as ^ ±00 


(2.28) 


where A and B are 0(1) constants. 

Outside the critical layer, the streamfunction consists of a WKB-type solution and 
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simple exponentials which satisfy Laplace’s equation, 


yc<y- 

i’+ 

-t 

1 

d 

II 

r~M^r(y-yc)+\iKy-ycM\h 

(2.29) 

kM{y-ycY 

0 < 2/ < 2/c : 

0_ 

= C-e-^y + D_e'^y+ 

r~M^r(y-yc)+\iKy-ycM\h 

(2.30) 

kM{y-ycY 

where C ±, D- 

and E± are constants whose order of magnitude will be set by matching to 


the critical layer streamfunction. The fact that '0+ and ^0+ /dy must vanish as y —>■ oo is 
imposed by including only the decaying exponential for y > yc- the WKB mode vanishes 
as 2 / —>■ oo when rj > 0 provided dr > 0. 

The governing equation (12.201) is third order and so 'ip,'ipy and 'ipyy must be continuous 
everywhere. The oscillatory form of 0^^ in the critical layer means it must match entirely 
to the outer WKB-type solution 

^WKB d _ ^-[Sdr{y-yc)+2iKy-yc) ]/ri ^2 31 ) 

k‘^{y-ycY 

either side of the critical layer so E± = This means that the outer WKB solution 

only contributes at 0 {y/fj) to the tangential velocity u near the critical layer whereas the 
critical layer streamfunction 0 forces an 0(1) tangential flow (see (12.271) 0 As a result 
there must be a large-scale flow of 0 ( 1 ) in both u and v outside the critical layer: this 
explains the scaling of the integration constants in (12.281) and means that C± and 0_ 
are all 0(1). There are then six (complex) conditions to be satisfied at leading order 
by the five remaining constants and complex eigenvalue a. The first four are matching 
conditions on u ( 0 y) and v (—ik-ip) as y —>■ 2 / 0 , 


- kc^e-^y'^ = 


— -f A & C+e-^y’^ = B 
2 ik 


(2.32) 


and y ^y~ 


- kC-e-^y’^ + kD_e'^y‘^ =-J—+A & C-e-^y’^ + D.e^y'^ = B (2.33) 

2 ik 


or eliminating A and simply two jump conditions across the critical layer 

/ 27r 


= 


& [?^]_ = 0. 


(2.34) 


The two remaining (boundary) conditions, v{x,0,t) = 0 and u{x^0,t) = 0, require 


C_+D_=0 & - kC- + kD_ - My = 0 

kyc 


(2.35) 


where, while the outer WKB-type solution contributes at 0{^) to u near the critical 
layer, it must contribute at 0(1) to u at the inflow boundary. The resulting dispersion 
relation is 


Say„-\ikyl ]/rt /^g-fcyc _ Q 

kyc V ik 


(2.36) 


which leads to the same leading expressions for dr and < 7 ^ given in (12.121) and (12.141) . The 
necessary change in the scaling of the WKB solution contribution gives the dominant 
0 (r/logl/r/) contribution to dr and the exact numerical counterbalancing of the large 
scale tangential flow gives the subdominant 0 (? 7 ) contribution. 















11 


Instability Driven by boundary inflow 


2.1.3. The Need for Shear and Inflow 


To emphasize that shear is a key ingredient of the instability, the above problem can 
be solved for the shearless flow U = x + Tyy leading to the requirement that there must 
exist 5ie(tT) > 0 with 


^-ky-ay/r,^- ^ q 


(2.37) 


which is never satisfied. The presence of an inflow boundary is also crucial: converting 
the above inflow boundary to an outflow boundary {rj —>• — 77 ) removes the instability. 
This is why the instability discussed here is not relevant to the considera ble literature o n 
suction boundary layers where suction is always a stabilizing effect (e.g. Joslin ( 1998h b 


2.1.4. The Third Boundary Condition 

Imposing the parametrised boundary condition (1 — fd)u + fidu/dy = 0 as the second 
boundary condition at y = 0 gives the modified dispersion relation 



(2.38) 


(/3 = 0 recovers non-slip and /3 = 1 a no-vorticity or equivalently stress-free condition 
as u = 0 along y = 0). For the boundary layer instabilities, the LHS is 0{^rj) so this 
dispersion relation can only be assumed similar to the non-slip relation (I2.15P if /3 <C y/rj- 
In fact, numerical computations indicate that the boundary layer instability is suppressed 
by /3 « 2.Qy/ri/k -C 1. For the critical layer instabilities, the condition (12.1011 becomes a 
possible balance between 3 different terms 


1 

a 



/{2kri) 


P 

1-/3 


(2.39) 


where cr now has an 0(1) frequency as in (I2.11|l . If /3 ^ 77 , the dominant balance is 
between the 2nd and 3rd terms (as opposed to the 1st and 2nd for non-slip) and now leads 
to damped eigenvalues. It is then clear that for instability to occur, there should be no 
restriction on the normal derivative of the tangential velocity at the inflow boundary. In 
practice this means that the instability only really occurs for non-slip boundary conditions 
when 0 < 77 <?; I which incidentally is the one choice which, in concert with the no-normal 
flow condition, means no disturbance kinetic energy is being advected into the domain 
through the inflow boundary. 


2.2. Half Plane: Viscous 

The base state (12.11) is unchanged (by design) when viscosity is introduced but the lin¬ 
earised disturbance equation ( 12.211 now includes diffusion of vorticity: 






Re 




(2.40) 


Looking for a modal solution w(x, y, t) = u}[y) exp( 7 fcx -I- at) leads to the equation 


d'^Cj 

dy^ 


Re Tj 


dw 

dy 


[Re{a + iky) -|- fc^)] w = 0 


(2.41) 


which has the bounded solution (as y —> 00 ) 


UJ 


el^^^VAi 


Rea + kf + jRe^rj'^ 
fikRe)’^/^ 


+ {ikRe)^^^y 


(2.42) 
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where Ai{z) is the Airy function bounded as \z\ 
relation is then 


oo with I arg( 2 )| < tt. The dispersion 


f 




Rea + + jRe^r]'^ 

{ikRey/^ 


(ikRe)^^^y 


dy = 0 . 


(2.43) 


and instability (5Re((T) crosses through 0 to become positive ) o ccurs at a critical i nflow 
Vcrit{k, Re). This integral is actually the same as that treated bv iGallet et al. ( 20101 1 (see 
their expression (39)) after the transformations 


/u2/3 

« ■ (fcAe)i/3 ’ 


4i?e 


N„k^ J 


1 


1/3 


{a + krf) 


(2.44) 


{Ng and Og from Gallet et al.l ( 2010h l. For the k of interest (< 0(1)), ky -C \a\ and we 


can reuse their critical value of Ng which has Og passing through the imaginary axis to 
give 


Vcrit = {\Ngk 


1/3 


Re 


-2/3 


^crit — 


4i?e 


1/3 




(2.45) 


where Ng = 4.57557 and Og = 5.63551/ (consistent with Gallet et al. ( 2010h who quote 
‘4.58’ and ‘5.62i’ for N„ 


‘‘9 

and a„ 


respectively). This threshold rjcrit tends to zero as k 


0 albeit with the unstable eigenfunction extending a distance 0 {{kRe)~^^^) in the y 
direction. In terms of connecting this analysis to other problems, there are two notable 
cases: k = 0 ( 1 ) which is the interesting case in rotating flow where the wavenumber is 
forced to be an integer by periodicity, and k = 0 {Re~^) which is gives the most uns table 
disturbance in a domain bounded in y (i.e. a channel see lNicoud &: Angilella (1199711 1. In 
the former case, rjcrit = 0 {Re~^^^) and the implication from the scaling of the critical 
frequency is the growth rate away from criticality (i.e. jiy — ricrit\ = 0 {rjcrit)) will be 
0(i?e“^/^) or 0 ( 1 / 77 ) as befor e. For k = 0{Re~^), rirrit = 0{Re ^) which is consistent 
with the numerical findings in Nicoud Sz Angilella I (|l997ll that the threshold ‘crossflow’ 
Reynolds number for inflow instability in their plane Couette flow is independent of the 
shear Reynolds number. 

For any given fc, further boundary-layer-type instabilities exist as 77 increases with the 
first 6 thresholds listed in Table 2 for fc = 1. Within this ‘boundary-layer’ scaling of 
A := —Re^^^Oi and f\crit '■= Re^^^ycrit for k = 0(1), asymptotic predictions for A —>■ 00 
can be derived following the same route as in the inviscid case. This proves a little more 
involved leading to two coupled relations 


Vcrzt = ^( Y log(27rA/A:) - - log 


A2 


Vcrit 


= ^Tr{8n — 1 ) 


.3fc log(27rA/A:)_ 
77 = 1, 2,3,... 


-1/3 


(2.46) 

(2.47) 


which work reasonably well for A = 0(10) given that higher order terms may only be 
log(A) smaller. 

As in the inviscid case, there are also interior critical layer modes excited for even 
higher inflows: if the critical layer is at y = —ai/k = 0 ( 1 ), then 


^crit — \^i\ 




(2.48) 
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asymptotic 


actual 


n 

A 

rjcrit 

A 

rjcrit 

6 

35.2992 

16.8776 

35.183 

16.501 

5 

29.5968 

14.2990 

29.482 

13.885 

4 

23.8395 

11.6711 

23.730 

11.201 

3 

18.0097 

8.9777 

17.9033 

8.4197 

2 

12.0754 

6.1886 

11.9737 

5.4762 

1 

5.9596 

3.2301 

5.8938 

2.1875 


Table 2. The 6 lowest inflow thresholds for instability in the dispersion relation (12.431) with 
fc = 1 as Re — >■ oo and asymptotic estimates from (12.461) & (12.471) . 


2.3. Inflow and Suction Together: Inviscid Plane Couette flow with Suction 
The half plane system exhibits boundary inflow instability for all 77 > 0. This can be seen 
by a simple rescaling of space 

K = rf'^^k, Y = y/rf'^^ & e = \l'rf (2.49) 

where /3 > 0 so that the equation ( 12.2011 becomes 

(cr + irSY + edY)iipYY — = 0 (2.50) 

which is just the original equation with a new small number e as 77 —> 00 . However, this 
is rather artificial since the applied pressure gradient also has to be increased with 77 to 
maintain the constant shear in y. Resorting to a constant pressure gradient instead now 
means the shear field decreases as 77 increases again making it difficult to draw general 
conclusions for large 77. Introducing another boundary is then the only alternative and 
this must be an outflow boundary if the resulting base flow is to be steady and spatially 
non-developing. The simplest modification is to add an outflow boundary at 7/ = 1 
which is moving at x so t hat the constant-vorticity basic flow m is still a solution 
( Nicoud fc Angilella 1997ll . The equivalent expression to (12.81) is then 


ifkiy) = - 


sinh(fc?/) 


sinhfc(l-7/)e-('"^+2*'=!'')/''dy 

sinh(fc7/)e~^'^^~'’2®^^ ^^'^dy 


k sinh k 

sinhfc(l — y) 


'V 


k sinh k 


with the dispersion relation 


/ sinh k{l - y)e-(‘"^+= 0 . 

Jo 


(2.51) 


(2.52) 


when non-slip is applied at the influx boundary y = 0 for rj > 0. This is essentially the 
same as the half plane dispersion relation and will have unstable eigenvalues as there is an 
inflow boundary. The simple rescaling ( 2.49L however, is disallowe d and instability is lost 
if r] becomes too large (see figure 12 of 


Nicoud fc Angilella ( 1997ll l. This could be inter¬ 


preted as the stabilising influence of the newly-introduced suction boundary ultimately 
overpowering the destabilising inflow boundary. Further evidence for this comes from the 
pressure-gra dient-free version of this flow which is also linearly unstable in the presence 
of viscosity ( Doering et al. llioOflll . In this case, the base state varies exponentially in the 
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cross-stream direction (see eqn 2.13 of lPoering et al. ( 200n[)i and possess es an area of 


linear instability in the {tan 0, Re) plane (figure 3 of lPoering et al. 


(20o 3) where tan 9 

is their proxy for rf). The lower boundary of this instability region, 6 iower{Re)^ plausi¬ 
bly scales like Re~^ which is the viscous threshold for the instability as described in >12.21 
whereas the upper boundary h as 9„.nrtpr(Re) = cot“^ 54,370, which appears to be suction 


ultimately stabilizing the flow (iHocking 1119741) . 


3. 2D Swirling Flow with Radial Inflow 

We now add curvature to the discussion and consider the basic, steady, axisymmetric 
solution 

V = --r + rn{r)e (3.1) 

r 

between two boundaries at r = 1 and r = a (> 1) with P(l) = 1 (i.e. the inner radius r* 
and the angular velocity there ft* are used as length and inverse timescales respectively) 
and 0 < ry so that there is a radial inflow. The 2P Navier-Stokes equations for the 
deviation u of the total flow Utot from the basic solution dSU, 

u := Utot - U = u{r, 6 , t)r + v{r, 6 , t)9 (3.2) 


are 


d ^ d rj d 

-h ft - - — 

dt do r dr 


r]u 


— 2flv -I- u. Vu — 


dp 


r dr 
1 


Re 


- 


r2 Q0 


d ^ d rj d \ dirfl) vv ^ „ uv 1 dp 

ft— - In -b - -V + f^u+u.Vv + — + -7^ 

r^ r r dO 


dt d9 r dr J dr 


1 

Re 




together with incompressibility 


where 


1 d{ru) Idv 

-i —- -I-=0 

r dr r do 


Re := 


r*'^ft* 


2 du 


and V is the kinematic viscosity. Introducing a streamfunction 

1 dij) dip 

reduces the system (lT3ll-(l^ to 

^ \ ^ V ^ \ \ ! I dZ dip 1.2/ 1t//a/\ 


where the Jacobian is defined as 

J{A,B) 


dA dB dA dB 
dr do do dr 


(3.3) 

(3.4) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 

(3.9) 
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Z := 


1 d{r'^ri) 
r dr 


(3.10) 


is the vorticity of the basic flow (EU. There are two special cases where the gradient of 
the vorticity vanishes: uniform rotat ion (Q = 1 so Z = 2) an d irrotational flow 17 = 1/r^ 
with Z = 0 originally considered by lllin fc Morgulis I ()2013ll . 


3.0.1. Basic state 


There are two usual scenarios: a) the swirl field 17(r) is set by the boundary conditions 
as in Taylor-Couette flow or b) the swirl field is determined by an imposed body (gravita¬ 
tional) force as in the astrophysical context. In the former, the swirl is determined by the 
azimuthal momentum equation given a radial flow with the radial momentum equation 
determining the pressure field. For example Ilin & Morgulis (2013) discuss the inviscid, 
irrotational, axisymmetric Taylor-Couette-like flow 


17 


inviscid 


1 

fa 


(3.11) 


which is the only possibility with axisymmetric radial flow (recall 17(1) = 1 by nondimen- 
sionalization). Gallet et al. (2010) consider the viscous Taylor-Couette situation where 
possible base flows form a 1 -parameter family 

VLtc = (3.12) 


with A a constant set by the motion of the outer cylinder and is generally rotational 
{Z = A{2 — riRe)r~'^^^). In the (latter) astrophysical context, the acting gravitational 
(body) force sets the swirl field (via the radial momentum equation) which then sets the 
radial flow by the need to balance ensuing azimuthal viscous stresses. The focus here is 
on the latter situation and we consider the general set of profiles 


17 = r“ 


(3.13) 


in order to understand how robust the boundary inflow instability is. Profiles with 
—2 < a < 0 hav e angular momentum increasing with radius and so are Rayleigh-stable 
( Rayleigh 1917ll . The gradient of the vorticity dZ/dr = Qf(a-|-2)17/r also does not change 
sign across the domain r € [!,«] so that the flow is inviscidly stable (for disturbances 
which va nish at r = 1 an d a) by a rotating flow analogue of Rayleigh’s inflexion-point 
theorem I Rayleigh I 88 OII . Particularly interesting choices for a are a = —3/2 which is 
the Keplerian profile for a thin accretion disk due to the radial force balance 


- rl7^ = - 


GM 


(3.14) 


(where G is the gravitational constant and M the mass of the central object) and a = —1 
which is used to model spiral galaxies. Since we work with deviations away from the 
basic state, the exact body force required to maintain the underlying rotation profile is 
not explicitly needed in what follows. In contrast to the Taylor-Couette situation, the 
azimuthal component of the Navier-Stokes equations forces the existence of a small radial 
flow 






a 1 

Re r 


(3.15) 


which is an inflow if dfl/dr < 0 (a < 0). Studying the consequences of this small accretion 
flow is the motivation for this work. 
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o-j 


Figure 2. Inviscid 2D instabilities for a = —2, m — 1 and a — 2 with scaled growth rate 
ar = (Jr/ (a) /(2a) plotted against frequency ai. The right (left) vertical dashed line is 
(Ti = —mQ{a) {ai = —mfl(l)). The eigenvalues from a full 2D eigenvalue calculation are shown 
as (the upper) black dots for rj — 10“® {N — 1000) and (the lower) blue dots for r] = 10“"^ 
(since N = 4000 here, the 15 most unstable eigenvalues are marked with a normal-sized dot 
and the rest by smaller blue dots). The red squares are the 20 most unstable modes from the 
boundary layer asymptotics (section 13.111 with the dash-dot red line tracing the path of the rest 
(eigenvalues calculated from the boundary layer equation (12.2211 are indistinquishable from the 
asymptotics at this scale). The green dashed lines are the critical layer asymptotic expression 
(13.2611 with r] — 10“® and 10~^ plugged in and appropriately rescaled: dr = 5iar,J— ^ and 

Gi = —mfl{Rs) with Rs taken over [l,a]. Notice that at 77 = 10“"^ even N — 4000 struggles to 
fully resolve the critical layer eigenvalues near the inner radius - see the hump in the numerical 
data which breaks the otherwise excellent correspondence with the asymptotic prediction. This 
hump is delayed to lower ai if N is increased until it eventually disappears - e.g. the rj = 10“® 
curve. 


3.1. 2D Inviscid Swirling Flow with dZ/dr = 0 

We study the simplest case of vanishing vorticity gradient in the base flow first {dZ/dr = 
0), so n = l/r^ or 1, to show how the analysis mirrors that in the half plane case. The 
case of uniform rotation fl = 1 initially looks uninteresting because the base flow needs 
an azimuthal as well as radial body force to maintain it (since a = 0 in (13.151) ). But it is 
worth considering as then only the enforced radial inflow creates shear in the base flow 
and the question is whether this is enough to generate instability. The inviscid, linearised 















real(u),real(v) real(u), real(v) 
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0.5 
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-0.5 


1 1.2 1.4 1.6 1.8 2 

r 




Figure 3. The most unstable boundary layer eigenfunction (upper plot) and an unstable critical 
layer eigenfunction with at ~ — 0.5 (lower plot) for rj = 10 “^, m = 1 , a = 2 and a = —2. In 
both, the green dashed line is Ke(u) and the blue solid line is Ke(n). The critical layer position 
is shown as a red dashed line for the critical layer eigenfunction. 


governing equation (13.81) is just 


d ^ ^ y 9 \ n 


(3.16) 


which is the ‘curved’ analogue of (I2.2|l and again 3rd order rather than the usual 2nd 
order. As before, the solution can be written down using characteristics as 

uj{r,9,t) = uJo( r^ + 2r]t,9 + — [ _ D,(r*)r*dr*] (3-17) 

\ V Jy/r^+2rjt ) 


where a;o(r, 0) is the initial vorticity. After a discrete Fourier transform, the e™® 
ponent is 


com- 
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where g is an arbitrary function. For modal growth, ipm (r, t) should only depend on t 
through an exp(tTt) factor for some complex constant a which forces 

giO = Bexpf^[ n{r*)r*dr*\ (3.19) 

V2?7 V Jvi J 

{B some constant which can be set to 1 providing the influx boundary conditions don’t 
force B = Q) and letting t) = exp (at) then 


1 i9 d m? 
r dr dr r^ 


This has the solution 



/ err^ 

= exp 


' a f'^ 

im 


V J 

, ( 

af^ 

dexpl 



n{r*)r*dr* 


i’mir) = G+{r,r;a,m)exp^^^ - ^ J n{r*)r*dr*^dr 

+ [ G-{r,f]a,m)exp( — -— [ fl{r*)r*dr*'\df 

Jr \‘2v V Jf ) 


where 


G-{r, r; a, m) 
G+(r, f;a,m) 



(3.20) 


(3.21) 


make up the Green’s function built around the 2 (usual) no-normal-velocity boundary 
conditions that t/’m(l) = 0 and ^prnio.) = 0- Imposing the third boundary condition of 
non-slip, d'4>m(fl)/dr = 0, gives the dispersion relation 




,, , err im 


VL{r*)r*dr* ) dr = 0. 


(3.22) 


For Q. = 1/r^ this is exactly expression (5.3) from lllin fc Morgulis ( 2013 1 which they 
establish has unstable eigenvalues. The energy budget for an infinitesimal 2D disturbance. 


d 

dt 






r—1 


du — ( ruv—— 
dr 


(3.23) 


where ( ) := rdrdd, u = 0 at r = 1 and u = v = 0atr = a, clearly shows that 

enhanced growth rates of 0 {y/rj) are only possible if dQ/dr is non-zero somewhere in the 
domain, emphasizing the importance of azimuthal shear. In particular, for D = 1, the last 
term drops so that any instability can only have a growth rate of 0{r]) at best. In fact, 
the dispersion relation appears only to have stable solutions, indicating that boundary 
inflow and just the shear due to the radial inflow are insufficient to drive any instability. 

The r; —>■ 0 asymptotic analysis for D = l/r^ mirrors that in the half-plane case with 
the saddle point at ct -I- imQ{rs) = 0. For the boundary layer scenario, adopting the 
scalings 


a = — imn{a) -|- ^/ryxja {dr + idi) + 0{ri) 


z := 




ja-r) 

Vd 


where x ■= —^rnil (a) (D (r) < 0 for cases of interest) and z is a boundary layer variable 
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(as before), retrieves (12.1511 which has the unstable eigenvalues a listed in Table 1: there 

is instability with growth rates 0 {(^ rj^—^mO,' {a)) / . 

For the critical layer scenario, the appropriate eigenvalue scaling is again (12.lip which 
given (T + imH{rs) = 0, leads to the expansion 

for the saddle point position where di = —mfl{Rs) defines Rs- The balance between the 
leading contribution from the saddle point at r = Ts (first term) and the endpoint r = a 
(second term) is then 


V—*?77r(r™+^ _ j,i m) / 


im 


-\mrsD.' {rs) 


exp - 

V 2?7 77 


I}{r )r dr ) — - -;-FTTr^expl 


{di + mfl{a)) 


da‘ 


= 0 


which leads to the frequency condition 

di{Rl - a^) 
2 r, 


Arg-( exp 
and growth rate 

Sdj- — 


-!!!:[ fi{r*)r*dr*+ - 

V Jr, 4 


= 0 


V 


(a2 - 


logfl) +21ogf(^^"-^""^)(7---—+0(1) 

\vj V {a^-a--)J-mRsn'{Rs) ) ^ . 


(3.25) 


(3.26) 


(a™ — a~^)^/^mRsD.' {Rs) 

Figure [2] shows that this performs well for numerical calculations with rj = 10“^ and 
10“^. The most unstable boundary layer eigenfunction for 77 = 10“^ is shown in figure [3] 
along with a critical layer mode with di ~ —0.5 which is typical. 


3.2. 2D Inviscid Swirling Flow with dZjdr 0 

The asymptotic analysis can easily be extended to treat more general rotation profiles 
fl{r) where the gradient of vorticity is non-zero. The idea here is to assume some small 
viscosity to induce a radial inflow (via[3T5|) and then to show that the inviscid instability 
mechanism is robust to the exact form of the azimuthal flow (e.g. whether it has a vorticity 
gradient or not). So, taking the linear, inviscid version of perturbation equation (13.Sp and 
inserting the usual ansatz 'if{r, 0, t) = ip{r) exp( 7777,0 + dt) gives 



f , . ^ Vd,\ im-ip dZ 

a H- %m\l(r) -L'?/’ =- 

\ r dr J r dr 

(3.27) 

where 

df Id m? 

(3.28) 


dr‘2- T dr 

The 77 = 0 problem 

(, + ,n.a(r))ct=''"*'‘f 

V / r dr 

(3.29) 


subject to the 2 boundary conditions 7/^(1) = '4’{a) = 0 seems to have no discrete modal 
solutions for the profiles = r°‘ studi ed here. This is easily proven for the special 
choices a = 0 & — 2 {so dZ/dr = 0) ( Ilin fc Morgulis I I 2 OI 3 II but is only suggested 
by numerical evidence generally. This observation is significant because it forces a non¬ 
standard singular perturbation analysis in which the additional flow component for 77 7 ^ 0 
has to contribute at leading order to help satisfy one of these 2 boundary conditions rather 
than just ‘hxing up’ the third boundary condition. 



























20 


R. R. Kerswell 


The analysis of the boundary layer instability is exactly the same as the SI = 1 jr^ case 
because the presence of a vorticity gradient does not effect the boundary layer problem 

at leading order. So the growth rate 5fte(tT) scales like {a)r]/a with azimuthal 

wavenumber m, local shear at the boundary — (a) > 0, and radial flow rj. The key point 
is that the boundary layer instability only depends on the shear at the inflow boundary. 

The analysis for the critical layer instability also proceeds in a similar fashion with one 
key difference: the problem for the large-scale flow is now complicated by a non-vanishing 
vorticity gradient but this only has a secondary effect: dZ/dr affects the growth rate at 
0 ( 77 ) rather than at the leading 0(77 log I/ 77 ) level. The starting point is the ansatz 

<7 = iai + Sdry/x/Rs (3.30) 

where = —mil^Rs) and Rg is the position of the critical layer, x ■= (Rg) and 

6 -C In the critical layer 

{-Sdr/y/rj + 2i^ + = 0 (3.31) 

to leading order where ^ := ■\yR^xJ^{r — Rg). This has the arbitrarily-normalised solution 

(3.32) 

Matching this to iprr outside the critical layer requires the WKB solution 

277 

(3.33) 

to exist there. There are two other ‘outer’ large-scale solutions either side of critical layer 
which solve 

[nir) - niRg)]C^ = (3.34) 

r dr 

and together accommodate the jump conditions 

M^ = 0 & [v]y = J^^^ (3.35) 

y ir] 

across the critical layer (found by integrating (I3.32p twice) and the no-normal velocity 
boundary conditions at r = 1 and r = a. Then, as in section [2.1.21 balancing the large- 
scale component of the azimuthal velocity with that from the WKB solution at r = a 
furnishes the growth rate and dispersion relation for the frequency. Only the higher 
order 0 ( 77 ) part of the growth rate depends on dZjdr or indeed 0 (r) whereas the leading 
0(77 log 1 / 77 ) part does not (note x which contains O {Rg) needs to be non-zero but 
otherwise scales out). 

3.3. Viscous Asymptotics for 0 < 77 <g; 1 for the boundary layer instability 

The inviscid asymptotics can be generalised to include viscosity which we do here just 
for the more dangerous boundary layer instability as this defines the viscous threshold 
for instability. In the absence of any other physics, the radial inflow is set by the viscosity 
present via (13.151) . However, to extract the scaling law for the viscous instability threshold 
for more generally-maintained radial inflows, we ignore this connection and assume 77 
can vary independently of Re. It is a simple matter to reinstate this connection later to 
establish stability or instability for a purely hydrodynamic situation. On this basis, the 




WKB 


-RsXd 


m?r‘^[Q{r) — ^(i?;,)]^ 


exp 


im 
V ■ 
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linearised disturbance equation for ip = is 


(T + imH{r) — 


ri d 
r dr 


imdZ 

Lip = — —ip- 
r dr 


Re 


C^ip 


(3.36) 


Setting cr + imLl{a) = etr with a = 0(1) and where e is the widt h of the boundary layer 
at r = a, and balancing the frequency, inflow and viscous terms ( Gallet et al. 201011 


t/e ~ 77 /e 


e'^Re 


(3.37) 


leads to e = Re and 77 = fje^ where 77 = 0(1). Defining the boundary layer variable 

(3.38) 


■v/X f a- 

with, recall, x ■= (a) and adopting the standard boundary layer decomposition 

ip = ipo + ipo + e{'ipi + ipi) + ..., leads to the boundary layer equation at r = a (the 
boundary layer at r = 1 is very weak and can be ignored at leading order) 




d^ipo 

dy^ 


= 0 


with 


N := 




X 


1/2 


& 


(To := a/{Nx^y^^- 


(3.39) 


(3.40) 


This is exactly equation (36) in Gallet et al. ( 2010ll if their ‘a’:= —(Tq. The solution which 
vanishes for 77 00 is 


dy^ 


iVl/3 




iV4/3 


(2i)2/3 4(27)2/3 


{2iY/^N^/^y 


(3.41) 


where Ai is the Airy function. Imposing the further non-slip condition dipo/dy = 0 at 
y = 0 defines the dispersion relation 


/ 


e2 




_ (2i)2/3 ^ 4(2i)2/3 


{2i)^/^N^/^y 


dy = 0. 


(3.42) 


The onset of instability is found for N = Nc where 5Re((7o) passes through zero (in 
the positive direction). Numerically, it is better to solve (13.391) directly as an eigenvalue 
problem rather than treat this integral equatio n (see Appendix A) . We Hnd Nc = 4.57557 
and (To = —5.635517 which is consistent with Gallet et al. (201^: the threshold radial 
flow rate for instability is 

Vcr^t = a/V//3yi/3^g-2/3^ (3 43) 

Since the viscously-induced radial flow in an accretion disk is 0{Re~^) <C rjcrit, the 
instability is not expected to be triggered unless it is substantially subcritical. Whether 
this is the case will be considered below in (JHl after we examine the effect of introducing 
some extra physics to the instability. 


4. 2D Compressible Swirling Flow with Radial Inflow 

From the astrophysical perspective, an important ingredient so far missing from the 
models considered is compressibility. Adding this extra physics also provides an oppor¬ 
tunity to test the robust of the boundary inflow instability. As the simplest model to 
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include compressibility, we consider the fluid to be an isothermal perfect gas so that 
pressure p and density p are simply related by p = c^p where is the isothermal speed 
of sound. The divergence-free flow 


U = -"^r+^e 


is a steady solution of the inviscid 2D momentum equations 

U— - — + — 

dr r po dr r^ 

+ ^ =0 
dr r 


(4.1) 


(4.2) 

(4.3) 


where poir) = p{r)/p(l) so the density is non-dimensionalised by its value at the inner 
radius r* and p^ir) = p(r)/(p(l)r*^r2*^). This non-dimensionalisation means po = 
where <5 := c^/(r*^D*^) and <5 —>■ oo is the incompressible limit. In a thin Keple- 
rian disk, however, the basi c rotation speed is highly supersonic so <5 <C 1 (e.g. p87, 
iFrank. King &: Raine (198^) and p ~ Re~^ -C S so the question is whether the insta¬ 
bility still operates for 77 <C d 1. 

In the laboratory, the body force on the right hand side of (14.21) would be absent 
leaving only the pressure gradient to drive the centripetal acceleration. In an accretion 
disk, however, the gravitational attraction to the central object plays this role and the 
pressure gradient only exists to balance the much smaller radial advection term. To model 
this latter situation, a body force is included (albeit here with different radial dependence 
so D := 1/r^) to support the irrotational basic state needed to satisfy (14.31) in the absence 
of viscosity. The pressure gradient is then only 0{p'^). With this, the density drifts very 
slowly as mass balance requires 


dpo 

dt 


y dpo 
r dr 



(4.4) 


Given the instability being studied here develops over a much shorter 0{\/y/rf) timescale, 
both this secular change and the very small pressure gradient can be ignored i.e. it is 
valid to set po = 1 and po = 43- After doing this, the 2D linearised disturbance equations 
are then 


^ ^ ^ ^ , on 

cr -I- im\L[r) -— Ml = — 5 -mi -|- Zilvi -— 

\ r dr J dr 

f , ■ nf \ ^ -7 

\(T + imil[r) - —\vi= -ttVi — Zui -pi 

\ r dr J r^ r 

{ Tj d '\ 

a + imQ{r) -— ]pi = —SV ■ Ui 

\ r dr J 


(4.5) 

(4.6) 

(4.7) 


where the disturbance fields Ui & pi are taken proportional to (pi 

Assuming the boundary layer scalings of 113.11 


cr = —im£l{a) + 


—mpil' (a) 
2 a 


^ = (a - J")! 


' —maVl' {a) 


2 p 


5pi). 


(4.8) 


{ui,vi,pi) 



—mail' (a) 



(4.9) 


t Formally, this is because ui.Vpo -C U.Vpi and pijpo dpo/dr <C 1/po dpi/dr. 














Instability Driven by boundary inflow 


23 


s 

Am((T) 

lRe(6') 

00 

-4.7795 

0.9114 

100 

-4.7795 

0.9114 

1 

-4.7796 

0.9092 

10-1 

-4.7784 

0.9290 

10"^ 

-4.7793 

0.9081 

10-® 

-4.7748 

0.9052 

10-1 

-4.0979 

1.1839 


Table 3. The most unstable eigenvalue a for the eigenvalue problem (I4.5l) - (l4.7ll with 7 ; = 10 
and various degrees of compressibility. The presence of S only becomes significant when it is 
< 0{rj). 


The equations become to leading order in rj 


0 = 2nv + ^, 

— ^rriD. {a)(a + 2i^ + 1 = —Z{a)u —^—p, 


d + 2if^ -\- 


df 

^ 5 ^ du ^ imv^ 


df 


ar] 


df a / 


which can be reduced to the equation 


, d\dv 


2aZ{a) 
mil' (a) 


a + 2^^ + £jp. 


(4.10) 

(4.11) 

(4.12) 


(4.13) 


This suggests that the incompressible limit (5 —> 00 which recovers the boundary layer 
equation (12.2211 ) still holds provided rj 6 and this is confirmed by numerical computa¬ 
tions: for example, see Table 3. 


5. 3D Linear Instability 

We now study 3 dimensional disturbances. Adopting the ansatz u(r, 9, z, t) = u( 7 .)g*(»n^'+fcz)+o't^ 
the linearized inviscid governing equations are 


^(7 + imll{r) 
a + imll{r) 


^cr + imilir) 


r] d 
r dr 
77 d ' 
r dr ^ 
r] d 
r dr 


u + ku — 2Ilv + 


dp 

dr 


rj im 

V - + Zu -\ - p 


w 


+ ikp 


1 dlru) im 

- - -1- V -b ikw 

r dr r 


0 , 

0 , 

0 , 

0 . 


(5.1) 

(5.2) 

(5.3) 

(5.4) 


In contrast to the 2D situation, discrete neutral modes start to emerge for fc yf 0 in the 
absence of radial flow. For example, in the case of axisymmetric modes the equations 
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Figure 4. The leading growing v eigenfunction calculated using (15.2211 in Maple for m = 1, 
k = 0.3, a = —3/2 corresponding to a = 0.64021—4.59526i (real/imaginary part is red solid/blue 
dashed line). 


(|5.1[) - (I5.4I) boil down to the 2nd order ODE for 


cPu 1 du 

dj-2 j. dr 


2 fc= 


zn + k^ 



(5.5) 


There are two profiles D(r) which make this just Bessel’s equation: D = 1 (a = 0) 
and n = 1/r {a = —1). In the former, unifo rm rotation case, eigensolutions are ax- 
isymmetric Poincare modes ( Greenspan 19681) . In the latter case, the general solution 
is M = AJ^{ikr) + BY^{ikr) where v := y^l + 2kP' /a'^ with the dispersion relation (since 
u(l) = u(a) = 0 ) 


Ji,{ik)Y^{ika) = Y^{ik) J^{ika). (5.6) 

This has purely imaginary eigenvalues a = iX with < 2k^. The issue is of course 
whether the emergence of discrete modes affect the instability. The answer is no as we 
now demonstrate again focussing on the stronger boundary layer instability. 


5.1. 3D Extension of the Boundary Layer Instability 

Working with the primitive variables, the appropriate scalings to capture the boundary 
layer instability are 


tr = —imLl{a) + i 


—mrjLl' (a) 
2 a 


f ■= (a-r)i 


' —maLl' (a) 


2rj 


(5.7) 


^/mrju, V, w 




(5.8) 


{u,v,w,p) 
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Figure 5. Leading inviscid 3D instability for a = —1.5, m = 1, a = 2 with scaled growth 
rate dr = Or/\J —^rnTjO! {a)/a plotted against /r, the degree of 3 dimensionality, as computed 

using the full eigenvalue code - y = 10”^ uppermost (blue) line with triangles, y = 10~® second 
uppermost (red) line with) squares, y = 10~® third uppermost (green) line with diamonds 
and y = 10”^ lowest (cyan) line with filled circles - where N = 1000 proves sufficient (The 
markers on each line actually show the N = 500 stability results for y = 0.4, 0.45 and 0.5 to 
demonstrate convergence). The leading instability is also shown for the boundary layer problem 
(15.141 & (15.151) as a dashed black line {N = 4000). The boundary layer scaling works well - 
the full eigenvalue prediction smoothly converges to the boundary layer prediction as 77 —^ 0 - 
for y < 0.3 but beyond this the scaling is no longer correct. This plot demonstrates that 
increasing 3 dimensionality acts to suppress the 2D instability. 


in the boundary layer ^ = 0(1) where 


y := 


k 


m 


(5.9) 
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measures the degree of 3-dimensionality. With these, the equation set jSIli-dEl reduces 
to 


0 = 20ti -I- ian'(a)^, 



w = 0 


at leading order respectively for small y'rmj. Eliminating u and p then gives 

-l- 2z^ ■ 


d\ dv 


d \ dw 




i"fv — 2iw, 


(5.10) 

(5.11) 

(5.12) 

(5.13) 

(5.14) 

(5.15) 


which is the 3D generalisation of the single boundary layer equation in 2D (see (12.221) 
where 

\/2{a + 2) 


7 := 2ap- 


{-a) 


(5.16) 


This reduced system only involves m through p and is valid for m -C 1/77. So, for p held 
fixed, the growth rate scales with ^/m (see the rescaling in ()5.7I1 ) across the astrophysical 
regime of 1 <C m l/?7. 

In the particular case treated by lllin fc Morgulis I (l2nl3^ (a = —2 so 7 = 0) adding 
an axial wavenumber does nothing to the 2D instability provided pk stays small. In the 
general case of interest — 2 < a < 0 and p > 0 (clearly the stability problem is symmetric 
under the transformation {p,v,w) —>■ (—/i, z;, —re) so only p, > 0 needs be considered) a 
solution is available in terms of Kummer functions. The equations (15.141) and (15.151) can 
be combined to a single equation for v 


where 


{g^ + 2 ig - 7 ^) 7 ; = 0 




(5.17) 


(5.18) 


This can be factored as {Q + 4iri)(^ -|- 4ir2)t' = 0 where Ei := i(l — -^1 — 7^) and 
r2 := 7(1 + v^l — 7^). Under the change of variable z := ji{d + 2i^)^, each factor is just 
Kummer’s differential equation, 


/I ^dv ^ 

^372+( 2 -^)—= o 




di 


J = l,2 


(5.19) 


The solution for v (and therefore indirectly w) which decays exponentially for ^ ^ 00 
and (Tr > 0 is 

V = e^AU{\ - Ti, i, -z) + BU{\ - r2, ^ -z) (5.20) 

where U is the multivalued Kummer’s function (see §13.2.25, NISlI ( 2010h ') and A and 
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Figure 6. Inviscid 3D instabilities for a = —1.5, m = 1, a = 2 and 77 = 3 x 10 ^ with scaled 
growth rate dr = Orj\J—mrfi (a) j (fla) plotted against frequency ai for various values of y,. The 
right (left) vertical dashed line is Ui = —mQ{a) {<Ji = —mfl(l)). The eigenvalues are calculated 
from a full 3D eigenvalue calculation with N = 2400 for (in order downwards from the uppermost 
2D eigenvalues): y = 0 (blue dots); y = 0.35 (red dots); y = 0.4 (magenta dots); y — 0.43 (green 
dots); y — 0.45 (cyan dots); y = 0.5 (black dots); y = 0.6 (blue dots with squares); y = 0.75 
(red dots with diamonds) and y = 1 (magenta dots with circles). This plot demonstrates that 
increasing 3 dimensionality suppresses the 2D instability (in fact completely hy y = 0.75 here). 
Notice that numerical errors start to creep into the eigenvalues whose eigenfunctions have critical 
layers far from the inflow boundary (the developing corrugations in the curves) as y increases 
(not unexpected as the numerical truncation N is only 0(1/77)). 


B are constants. If these are not both to vanish when the further boundary conditions 
V = w = 0 are imposed at ^ = 0, then the product 

C/(i - Fi, I,-\id^)U{\ - Fa, i ,= 0 (5.21) 


which defines unstable {dr > 0) eigenfrequencies d. The eigenvalues for 77 > 0 which 
smoothly connect with those at /r = 0 are given by f7(i — Fi, i, —^id^) = 0 since 77 —>■ 0 
implies Fi —>■ 0 and convergence to the 2D equation Qv = 0. These eigenfunctions can 
be plotted using 


Gamma(2) 
Gamma(l —r 1} ^ 


W(i - Fi, i -z) - e^C/(i - Fi, i -^) 0 < ^ 


< 


w 


ye^U{\-Ti,\,-z) -l{dr + di)<£, 

(5.22) 

and w = 2Tiv/y (with z := ji{d + 2iffj^ and ‘Gamma’ indicating the Gamma func¬ 
tion ) which compensates for the branch cut along the negative real axis in U routinely 
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Figure 7. The growth rate cr^ verses radial flow r] for linear 2D perturbations with mo = 1 
on the ID basic state ini) for a = —2il2, a — 2 and Re = 10®. The blue dots indicate the 
eigenvalues using 200 radial modes to represent each component of the perturbation velocity 
field and the red circles are a sampling of results of using 400 to confirm convergence. The leading 
perturbation (top curve) becomes unstable - Ur > 0 - at about Re = 5 x 10® and remains the 
only instability at Re — 10^. 


imposed by packages such as Maple (M is the single-valued Kummer function: §13.1.2, 
Abramowitz &: Steeun ( 1964ll ): see Figure |H 


The equations (15.141) and (15.151) can also be directly treated numerically (see Appendix 
A). This boundary layer approach works well for the dominant instabilities as /r increases 
from zero showing how they are systematically suppressed until their growth rates are 
comparable with the interior critical layer modes: see figure [5] which shows this for the 
leading instability. At this point (which is ^ ^ 0.4 for a = 2, r; = 10“^, a = —1.5) the full 
3D eigenvalue (I5.1l) - (l5.4p must be solved (see Appendix A) which shows the complete 
suppression of the all the instabilities by /i = 1: see figure [H 

The general conclusion is that 3D disturbances are less unstable than 2D disturbances 
under boundary inflow. In fact since the boundary layer instability depends only on the 
shear at the boundary, this could have been anticipated by usi ng a Squires tra nformation 
to map a 3D disturbance to a more unstable 2D disturbance ( Squires 19331) . Normally, 
such a transformation fails in a rotating system due to curvature terms but here these are 
marginalised by the dispersion relation being only sensitive to the shear at the boundary. 


6. Nonlinearity 

In this section we examine the nonlinear aspects of the boundary inflow instability. The 
natural starting point is a weakly nonlinear analysis of the 2D boundary layer instability 
in the presence of viscosity so that there is a finite (radial flow) threshold for instability. 
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Figure 8 . Upper: kinetic energy per unit length in 2 of the 2D perturbation, E, verses the radial 
flow 77 with the solid blue line representing the mo = 1 solution branch and the dashed red line 
the mo = 3 branch for Re — 10"*, a = —3/2 and a = 2 (resolution M = 20 and N = 60). Lower: 
the 2D streamfunction for the mo = 1 solution branch plotted (left to right) at 77 = 9.64 x 10“®, 
3.67 X 10~^ and 5.28 x 10“^ (shown as dots on the solution curve). In each, ten contours are 
plotted (dark/red-to-light/white being -ve to +ve) and at 77 = 3.67 x 10“^, maxi/ = 0.0478 and 
min 7 / = —0.0662 (in all, the outer colour/shading indicates the zero contour). 


Then the leading question is whether the new 2D solution branch exists for 77 < rjcrit 
(a subcritical bifurcation) or for 77 > rjcrit (a supercritical bifurcation). This question is 
most straightforwardly posed with the growing instability not permitted to change the 
(azimuthally-averaged) radial flow i.e. this is the control parameter of the flow. 

6.1. Weakly Nonlinear Analysis of 2D Viscous Boundary Layer Instability 
The analysis revolves almost completely around the boundary layer and can be handled 
relatively straightforwardly: see Appendix B. The growing instability is found to be 
supercritical so that the branch of 2D solutions exists for 77 > rjcrit- Along this branch of 
solutions, the azimuthally-averaged radial pressure drop, Ap, across the domain, which is 
a more natural control parameter for a laboratory experiment, changes. Furthermore, in 
an accretion disk, maintaining constant angular momentum of the flow is a more natural 
constraint under which to study flow bifurcations. These, however, just represent different 
perspectives of viewing the same bifurcation and the ensuing 2D branch of solutions. 
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Figure 9. Left: surplus azimuthally-averaged radial pressure drop 5p (blue dashed line), dis¬ 
turbance angular momentum 51 (red solid line) and surplus of 2D rotational energy above ID 
rotational energy 5E (black dash-dot line) all plotted against rj for the 2D (ttiq = 1) solution 
branch of £gure[ 8 l Right: the disturbance angular velocity profile for the 2D (mo = 1) solution 
at 77 = 9.64 X 10~^ Re = 10^, a = 2 and a = —3/2. Notice this closely resembles the prediction 
of the nonlinear analysis given that the boundary layer thickness is relatively large ( 10 “'^'^^) and 
the boundary layer oscillation extends 6 boundary layer thicknesses outwards. 


X 10" 




Figure 10. Left: r^uvdd against r for the bifurcating eigenfunction at 77 = 9.645 x 10”^ 

(left dot on (upper) figure | 8 ]). Right: the mean angular velocity profile for the mo = 1 2D 
solution at 77 = 3.67 x 10~^ Re — 10"^, a — 2 and a = —3/2 (solid blue line) and the profile 
D = r“ corresponding to the ID basic state (red dashed). This plot clearly indicates that angular 
momentum has been transported outwards by the saturated instability. 


For example, solutions with constant Ap and varying rj are the same as those with 
constant 77 and varying Ap provided the rest of the boundary conditions are consistent 
(e.g. azimuthally-asymmetric radial flow components vanish at the boundaries in both 
and the azimuthally-asymmetric pressure distribution on each boundary is unconstrained 
in both). To make this clear and to be able to make statements away from the vicinity 
of the bifurcation point, we now compute the fully nonlinear 2D solution branch. 
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Figure 11. The azimuthally-averaged radial pressure drop Ap is plotted again the azimuthal- 
ly-averaged radial flow rj for the undisturbed ID basic state (blue solid line connecting points 
A,B and E; Apio ■= (1 — l/n) + |i7^(l ~ 1/®^)) a-nd th® 2D solutions (red solid line connecting 
points A,C and D) which emanate out of the bifurcation for mo = 1, a = 2, Re = 10'* and 
a = —3/2 (as shown in Figures jS]). The bifurcation point is point A and the plot shows that 
the bifurcation is supercritical for either rj or Ap being the control parameter. Note that if Ap 
is fixed, the bifurcation leads to states with lower radial flow (compare states C & D with E). 


6.2. Fully Nonlinear 2D Solutions 

To study (non-helical) 2D bifurcations off the ID steady state (13.11) . a Reynolds number 
of Re = 10* is chosen as a compromise, being hopefully large enough to be in the 
asymptotic regime but not too large that flows become too arduous to follow numerically. 
According to the asymptotics, ijcrit = ~ 0.0061 (0.0087) for mo = 1 

(mo = 3) whereas actually the thresholds pcrit = 0.009635 (0.0136) are found numerically 
at Re = 10* (a = —3/2 and a = 2). Figure [7] shows how the spectrum of the linear 
operator depends on rj for Re = 10®. Instability is first possible at i?e « 5 x 10® (not 
shown) with the second mode of instability appearing for Re > 10* and by Re = 10® 
there are 3 unstable modes. Further computations for Re = 10® (not shown) show that 
further modes of instability emerge (now 8) and the persistent feature that each unstable 
mode is only uns table for a fi n ite ra nge of p. The lower threshold p; (which is the focus 
in this work and iGallet et al.l ( 2010ll l must tend to 0 as Re —?► c» to be consistent with 
the inviscid analysis whereas the upper thresh o ld r]„ must tend to a finite limit to be 


consistent with the analysis of lllin fc Morgulis (|201511 . 


The bifurcation is a Hopf bifurcation and the oscillation can be made to look steady 
by going into a frame rotating at the phase speed eg := aijmo at the bifurcation point. 
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Subsequently, the 2D solution branch is traced out by using the representation 


u 

N-1 M 

Umn (t ) ( r„+2 (C) - T„ (^) ) 



V 

-EE 

Vmn{t){Tn+2{0 - Tn{0 ) 

^immoiS-cgt) ^ ^ 

(6.1) 

_ p . 

n=0 m=0 

Pmn{t)Tn{Cl 




where ^ := (2r — a — l)/(a — 1), mg indicates the rotational symmetry of the bifurcating 
eigenfunction and ce is the azimuthal phase speed which is found as part of solution 
process. Figure [5] shows the new branch of 2D solutions arising from the one unstable 
mode present at Re = 10^ for mo = 1 and mo = 3 found by a Newton-Raphson root¬ 
ing finding algorithm with truncations varying from {M,N) = (20,60) to (10,150) (so 
0(10'^) degrees of freedom). Both bifurcations at the lower threshold rji are supercritical 
consistent with the weakly nonlinear analysis of ^6.11 and the solution branches reconnect 
with the ID solution (13.11) at the upper threshold rj^: see Figure [5] 

Along the 2D solution branch, the ‘surplus’ azimuthally-averaged radial pressure drop. 


r - 7/2 

6p := {Ap) 2 D - (Ap)iD = / 2il{r)v H- dr, 

Ji r 


( 6 . 2 ) 


(where (•) := l/27r is just an azimuthal average), the surplus of 2D rotational 

energy compared to the ID solution. 


6E := 27r / r( flv + ) rdr ) 


(6.3) 


and the disturbance angular momentum 51 all are positive quantities: see Figure O The 
initial increase in the pressure drop and the total angular momentum indicates that 
if either were used as a control parameter instead of the azimuthally-averaged radial 
velocity, the bifurcation would remain supercritical. For example. Figure [TT] shows a 
plot of Ap against the azimuthally-averaged radial flow parameter 77 for the 2D solution 
branch already plotted in Figured The weakly nonlinear analysis in Appendix B takes 
rj as the control parameter and finds that the bifurcated 2D solution branch exists for 
V > Vcrit (marked as the point A in Figure fTTI) with Ap larger for the 2D solutions than 
the equivalent (same 77 ) ID solution (see points B and C in Figure fTTI). If Ap is the 
control parameter, the bifurcation is still supercritical as 2D solutions can only exist for 
Ap > Ap at A but now the bifurcated 2D solutions give rise to 2 smaller radial inflow 
solutions (points C or D in Figure fTTI) compared to the ID solution E. 

Another interesting issue for accretion disks is whether this instability acts to transfer 


angular momentum / outwards. Forming the integral ' 
equation 


dO gives the conservation 


where I = f“l(r)dr and 


I(r) := / r[rfl{r) + u] rd9 I -|- = 0, 


rd9 


J := 


1 2 d ( V 

ruv — r]v ——r — — 
Re dr \r 


(6.4) 


(6.5) 


is the associated radial flux of angular momentum (the flux vanishes for the ID basic 
state (13.11) 1. The first (Reynolds-stress) term in J is computable using the bifurcating 
eigenfunction alone and can be used to determine the effect of the instability on the 
angular momentum distribution (the other two terms, which involve nonlinear aspects of 
the instability, subsequently balance the first term to give a finite amplitude 2D steady 
state: figure [9] (right) actually shows that r'^uv ~ rjrv). Figure [TO] plots this term over 
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Figure 12. The 3D solution branches continued out of the six Hopf bifurcations (see the inset 
which magnifies the dashed box for detail) from the 2D branch (thick blue loop as in figure [8| 
found for mo = 1 at Re = 10“^, a = —3/2, k — 1 and a = 2. All the bifurcations are supercritical. 
The ordinate is the disturbance energy (2D and 3D) and the abscissa the radial flow y. The black 
dot indicates the particular 3D flow state shown in detail in hgure [TJl Typically a resolution 
of {N, M, L) = (40, 20, 4) was used to follow these solutions with curve segments only shown if 
they are robust under truncation changes. 


the radius for the mo = 1 state just after the bifurcation (shown in figure H] as the 
leftmost point at 77 = 9.645 x 10“^) and shows that it alternates in sign but is mostly 
positive indicating outwards angular momentum transport. The angular velocity of the 
2D state at its maximum amplitude (the middle mg = 1 point in the upper plot and the 
middle cross-section in the lower plot of figure [5]) provides more definitive evidence of 
this outward transport: see figure ITUl 


6.3. Fully Nonlinear 3D Solutions 

To briefly explore the possibility of 3D states, 3D bifurcations were sought off the 2D 
branch of mg = 1 solutions. Concentrating on 3D states with an axial wavenumber 
fc = 1, six Hopf bifurcations were found between the initial 2D bifurcation point y = 
V 2 D = 0.0096 and 77 = 0.0152. These were then traced using the fully nonlinear steady 
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representation 

^i{mmo0-+lkz-) ( 6 . 6 ) 

by moving into a frame moving with the phase speed Cz in z (initially Cz '■= Ui/k at 
the bifurcation point where ai is the instability frequency). Since 6* := 6 — cgt already 
incorporates the phase speed of the initial 2D instability, these secondary 3D states are 
steady in a rotating and translating frame. Figure [T^] shows that all the traced bifurcations 
are supercritical i.e. there is no evidence of solution branches reaching r] < r] 2 D ■ Solution 
branches are shown as far as they are reproducible using different truncation levels. The 
maximum realistic resolution was (TV, M, L) = (40, 20,4) giving 58,963 degrees of freedom 
since the branch continuation approach was a direct Newton-Raphson solver albeit with 
multithreaded linear algebra software. A typical 3D flow is shown in figure [T3] with the 
presence of vertical jets at the outflow boundary and the lack of any large scale structure 
particularly noteworthy. 


u 

V 

w 

P 


N-l M L 

-E E E 

n=0 m— — M l—O 


Ulmn{t){Tn+2{0 - Tn{0 ) 
Vlmn{t){Tn+2{C) ~ ^n(C) ) 
Wlmn{t){Tn+2{C) ~ ^n(C) ) 
Plmn{t)Tn[C) 


7 . Energetics 

The instability is inviscid in nature so we first consider the energetics of this simplest 
situation: note that in the absence of viscosity, a = —2 gives the only consistent solution. 
If Utot = U + u is the total velocity field, then the scalar product of this with the 
governing Euler equations gives 

= ~ j l“L(uto4-d5') - j) putofdS (7.1) 


With no disturbance u = 0, this amounts to 


putot-dS' = 'n J dO J dz 


p{a,9,z,t)-p{l,9,z,t) 


= 7177 ( 1 + 77 ^) ( 1 - 


(7.2) 


so a net pressure drop radially inwards across the domain drives the radial flow and 
works at a rate to replenish the net kinetic energy leaving the domain (the rightmost 
term using U = — 77 /rf + l/r0). 

The need for interior shear to fuel the instability is apparent by looking at the distur¬ 
bance energy balance. Taking the scalar product of the disturbance velocity u and the 
disturbance evolution equation, leads to 


d 

dt 



— (u • VU • u) — 77 




d9 


or explicitly 



dfl 


ruv—— 

dr 



d9 

r—1 


(7.3) 


(7.4) 


so the disturbance can only gain energy through the underlying shear field (the last term 
on the right hand side is the loss of kinetic energy through the outflow boundary). The 
second term on the rhs is the energy transfer from the swirl field and this has to be 
positive (i.e. the underlying swirl field supplies energy to the disturbance) for a growing 
disturbance to achieve growth rates ;g> 0 ( 77 ). 
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Figure 13. A typical 3D state (shown as a dot on figure I12|l calculated using 
{N,M,L) = (40,28,4) (82,003 d.o.f.) plotted at 2 : = 0 (top left), 2 = 7r/2 (top right), z = n 
(bottom left) and 2 = 37r/2 (axial wavelength = 27r). The axial speed is shown by 10 contours 
going from -0.025 to 0.035 (the colour/shading of the zero contour is shown on the outside of 
the flow domain). The arrows indicate the speed and direction of the flow in the {r,6) plane 
with the longest arrow representing a speed of 0.272. 


Adding viscosity complicates the (energetic) situation by introducing interior dissipa¬ 
tion and the further possibility of viscous stresses at either or both radial boundaries 
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inputting energy into the flow: explicitly 




^277 


dO / dz 


S, z, t) - S, z, t) 

+ Ptot{a, 0, z, t) - ptot(l, d, z, t) 


If 1 

— ® Utot ■ 2e • dS —r— ( 2e : e) 
ReJ i?e' ' 


(7.5) 


where e is the rate of strain tensor for Utot- However, non-slip boundary conditions and 
constant r] do at least fix the net advection of kinetic energy out of the domain. Then 
either an increased radial pressure drop and/or increased work done by viscous stresses 
at the walls must offset the greater dissipation of a 2D bifurcated state. Certainly the 
nonlinear computations of ^ suggest that the radial pressure drop does go up for the 
bifurcated solutions. For an accretion disk, however, it may be more realistic to assume 
that the radial pressure drop is fixed but the radial flow rj can change instead. In this 
case, rather paradoxically, figure [TT] indicates that the 2 D flow will adopt a smaller radial 
flow to accommodate the greater dissipation of the 2D state. This runs counter to the 
usual argument that a turbulent disk should set up, via an enhanced eddy viscosity, a 
more accreting flow travelling down the gravitational potential to power it. 


8. Discussion 


In t h is paper we have revisit e d the 2D instability discu ssed separately bylNicoud fc Angilella 


(1997), Doering et al. ( 200fll l. Callet et al. ( 20in[ l and Ilin fc Morgulis (j2013l l 


ous contexts to explore the interesting mathematical and physical issues surrounding 
it. A simple half-plane model indicates that the instability operates by the crossflow 
advecting vorticity introduced by the inflow boundary across the cross-stream shear. 
Imposing vanishing vorticity at the inflow boundary eliminates the instability. This half¬ 
plane model also makes it clear that curvature or rotation is not important other that to 
restrict the allowed streamwise wavenumbers (azimuthal periodicity prevents very long 
wavelengths which are the most unstable in the rectilinear situation). It also highlights 
the fact that only inflow is destabilising suggesting that in situations where there is 
both an inflow and an outflow, the inflow boundary could be expected to destabilise 
the flow initially as the crossflow is turned on before the outflow (or ‘suction’) bound¬ 
ary eventually stabilises the flow at a higher crossflow value. This is clearly seen in 
circumstances where the u n derlying shear fl o w is not modified by the crossflow (e.g. 
Nicoud fc Angilella ( 1997 1. Ilin fc Morgulis ( 20I3h and here in but seems gener ¬ 


ally true even if it does (e.g. Doering et al. ( 2000 ): Fransson and Alfredsson (l2003ll: 


Gallet et al.l ( 20I0l) : Guha fc Frigaard] ~( 20I0l) : Ilin fc Morgulis ( 2015 ): Deguchi et al 


( 2014h although note Hains ( 197ll) ). The identification of the instability with inflow 


also explains why this instability is relatively unstudied compared to outflow boun daries, 
long lauded as a reliable means to stabilise unidirectional flows (e.g Ijoslin I (|l998h ). 

One of the most interesting aspects of this ‘boundary inflow instability’ is that the 
growth rates scale as ^ for boundary layer modes and r] log I /77 for interior critical layer 
modes when the perturbing crossflow is only 0{r]). This means that the instability has 
to draw energy out of the underlying shear field. However, the 2D nonlinear solutions 
computed in the rotating situation do not show any despinning of the swirl field but 
instead indicate that the radial pressure gradient which drives the crossflow more than 
replenishes this energy by working on the flow. 
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The apparently delicate mathematical structure underlying the instability - the lack 
of any discrete normal modes for the inviscid, vanishing crossflow situation - suggests 
that additional physics may easily suppress it. However, the work bv iGallet et al. (201^ 
alread y indicated that it could survive the addition of viscosity (see also lllin fc Morgulisi 
( 2015h ') and here we have also explored differe nt (more Rayleigh - stable) rotation profiles 
plus the addition of 3-dimensionality (see also Ilin fc Morgulis I ( 2015blH and compress¬ 
ibility. It is true that except for the rotation profiles, none of these have actually enhanced 
the instability (e.g. 2D modes are more unstable than 3D modes) but neither have they 
immediately suppressed it either (e.g. the viscous threshold fo r the instability is still only 
a very small crossflow of 0{Re~^) in a rectilinear geometry ( Nicoud fc: Angilella Ill99ffl 
and 0(i?e“^/^) in the rotating situation; see 113.3|1 . The conclusion is therefore that the 
instability is relatively robust except to the exact boundary conditions imposed at the 
inflow boundary and one should therefore expect instability whenever there is boundary 
inflow together with shear. 

Weakly nonlinear analysis of the primary bifurcation and branch continuation of the 
fully nonlinear 2D solutions indicate that the instability is supercritical in the cross- 
flow. Furthermore, secondary bifurcations to 3D states also appear supercritical (6 were 
found and continued) suggesting a succession of bifurcations in which the flow gradually 
becomes more complicated spatially and temporally as the crossflow is increased. This 
all makes the boundary inflow instability an inviting target for an experimental study 
especially as the primary instability is oscillatory and hence clearly identifiable. 

Astrophysically, we have established that a Keplerian rotation profile with crossflow 
of —rj/rr and non-slip boundary conditions at the inflow boundary will be unstable if 
the crossflow rj > 0{Re~^^^) - hence Rayleigh’s stability criterion can be circumvented 
by crossflow. However, in a quiescent disk, (molecular) viscous stresses only generate 
a crossflow of 0{Re~^) so the linear instability is not triggered. Moreover, both the 
(2D) primary and (3D) secondary instabilities are supercritical so there is also no ap¬ 
parent opportunity to reach the instability via finite amplitude perturbations at such 
low crossflows. Even if the crossflow was large enough, one would have to argue that the 
appropriate boundary conditions were in place at the inflow (outer) boundary of the disk. 
Nevertheless, the analysis presented here does highlight the potential for mass entering 
a disk to disrupt the orbiting flow if this mass flux possesses vorticity and also showcases 


the co mplications of introducing inflow boundaries in disk models (e.g. iKersale et al. 

(12004 1. 
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Appendix A. Numerics 

The 3D boundary layer equations (15.141) and (15.1511 can be directly treated numerically 
by transforming the domain ^ € [0, oo) to a; £ [—1,1) via the definition a: := {f—L)/{f^+L) 
where L is a scale factor {L = 1 works well here), v and w are expanded as differences of 
consecutive even or odd Chebyshev polynomials 

N 

(u,u)) = ^(a„,&„)(T„+2(A)-r„(A)), 

n—0 


where T„(x) :=cos(ncos ^ x) (Al) 
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so that V and w are forced to vanish at ^ = 0 and ^ oo, and the two equations 
(j5.14|) and (j5.15|) are collocated across the N zeros of Tn{X). The resulting 2N x 2N 
eigenvalue problem (with N varying from 100 to 10,000) is then solved using LAPACK 
routine ZGGEV (solving the 2D boundary layer equation (I2.22p is just a N x N special 
case of the 3D problem). The viscous boundary layer equation p.39l) can be similarly 
solved. 

The full 3D eigenvalue (I5.1l) - ()5.4p can be solved by expanding the primitive variables 
as follows 

Vr„+i(F)-r„(y)),d„r„(y) 

Cn J 

(A2) 

where Y := (2r — a — l)/(a — 1) so that the appropriate boundary conditions - u = 
0|r.=i, u = V = w\r=a = 0 - are automatically imposed. This 4A^ x 4A^ eigenvalue 
problem is again solved by LAPAGK routine ZGGEV typically with N = 2400. An 
associated inverse iteration code was developed where N could reach 12, 000 on a 128GB 
machine to confirm eigenvalues. 


V 

w 




-E 


,(T„+2(y)-T„(F)), 


Appendix B. Weakly Nonlinear Analysis of 2D Viscons Bonndary 
Layer Instability 

We introduce two small quantities: <5 as the amplitude of the instability and e = 
as a measure of the viscosity. The analysis proceeds from the following expansion of the 
velocity field 


^ = '^1 (d™’ (’■) + ^lo'' (y) ) + e( V’lTl*') + dr V)) + • • • I *) 



+■5(1) (^)wy v/={(dS”’M +dS’"’(!/) ) + 

+ (^) ^d®x"^:^|(do‘’(e + d(rly)) + c(e)|£(^,;) 

+<(4)' (^) «yy-=«{(dr’(-')+dr’(»))+ o(.)}e(9, 0 “ 


+ ... + C.C. (Bl) 


Ei0,t) 


—imQ{a)-\-ed’ )t 
^ 5 


a := ^ d2 + ..., 

r^:= + + 


where 
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with Nc = 4.57557, ao = —5.63551i and fjc = being the critical values (see 

1 13.3[l and x = —^rnD, (a) (only flow components in red need to be calculated). The 
objective here is to look around the bifurcation point rj = e'^fjc where 5 = 0 (fixed Re) to 
find a value of 7)2 consistent with 5 > 0: 7)2 < 0(> 0) indicates a subcritical (supercritical) 
bifurcation. We work with the full 2D disturbance equation (13.81) 


9 X 9 V 9 \ . . 1 dZ dtp 1*2/ /T-.r,\ 

S + - 7 *)^^ = 717 w + ife* (B2) 


B.l. Linear problem at 0(6) 
The boundary layer problem is 4th order in y 




d'^tp 


(m) 

10 


d7/2 


= 0 


(B3) 


with the boundary conditions that tp^f dtp)f^’/dy and d'^tp'f^’/dy"^ all vanish at large y 
and dtp^'^/dy = 0 at 7 / = 0. The outer (leading order) problem is 2nd order in r 


(m) 


72T,(m) 


f2(r) — rt(a) 


^10 - ^ dr 


(B4) 


(m ^ 0) to which the appropriate boundary conditions are tp^\l) = 0 ( 77 ( 1 ) = 0) and 
(®) = (Formally there is a weak viscous layer at r = 1 of thickness 

0(l/i?e) where the streamfunction boundary correction is 0(l/i?e) to accommodate the 
non-slip condition 7 ;( 1 ) = 0 but this plays no role in the nonlinear equilibration of the 
instability centered on the other boundary and will be ignored here and below.) The 
boundary layer problem alone identifies the bifurcation point (see 


B.2. Nonlinearity 0(5(5/e^)) 

Nonlinearity comes in at 0(5(5/e^)) generating boundary flows tp^^ and 7 ^ 20 ™^ follows. 
For 75 ^ 0 ^ 


^ ^-N — 
dy^ “ dy^ 




^10 




*(?7l) 

10 




^10 




dy dy'^ 


-I- c.c. 


(B5) 


where the solution which vanishes as 7 / —7 00 is sought. In the boundary layer, the outer 
solution 7 / 7 ^™^ is just the constant —tp[^\o) and since iSte(i\dip[^'^/dy\'^) = 0, equation 
dEH) can be integrated twice to 


/V ^ 
57/2 c dy 


and the required boundary layer solution is 


= fo(y) ■= 25ie<{ -i(tP['^\y) -tp'i^’(0)) 


fo(x) dx. 


dy 


^2o\y) = 


i 


°° X _ gNciv-x) 


(B6) 


(B7) 


This cannot be made to satisfy the no-slip condition dtp^^ /dy = 0 at y = 0 and so drives 
an interior mean flow tp^}_^ via 


dtp^°}^ 


dr 


d'fw 


dy 


= 0 . 


(B8) 


y=0 
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to ensure u = 0 at r = a. The interior mean flow problem for 




is just 


dr \dr r J dr 




= 0 


with b.c.s 




dr 


= 0 , 


#2°-l 


dr 


(m) 


since the leading interior nonlinear term J('!/’io 
at an 0(e) smaller level. The solution is 


A ; *(” 

^i’lO 




7o(x)dx 


(B9) 

(BIO) 


) + c.c. only drives a mean flow 


#2°-l 

dr 



1 ) 

1 ) 



e ^‘=^fo{x)dx 


2539- 


1 ) 

1 ) 


(Bll) 


so the mean flow (oc —dtjj^}_i/dr) is a decreasing function of the radius until the boundary 
layer is reached. There the mean flow undergoes an oscillation finally increasing sharply 
close to the boundary: see figure [TTl 


For 

2 7 i'^) ■ 

+*(^ro(2/)-^ro(o))^^ 

(B12) 

A particular integral for /dy'^ can be generated by the variation of parameters 

method given that the complementary problem, 






d^ll). 


20 


d ?/2 


_d 

dy 


—i 


dy 


^ - 2a-^ - -a^ + j^y) 

dy^ dy 


dy'^ 


= 0 


(B13) 


has solutions + jy) and e“^i3z(/? + yj/) but the ensuing integral expressions 

become unwieldy when integrated twice to get '020™^ which is needed below. Instead, it 
is better to numerical solve (|B 121) directly. In fact (IB 121) can be integrated once to 


dy^ 


(f d 

- - 2iVc(d-o + 2i?/) — 

dy dy 


■AiN 




(2m) 

20 


= /2(y) := 


—i 


dj^ 

dy 




iM. 

dy^ 


(B14) 


where a vanishing solution is sought as j/ —^ oo. The (numerical) solution strategy is 
to actually work over the domain y G [0,ymax] (with ymax ‘large’ but finite) which is 
transformed to a: G [—1,1] by the definition x := 2y/ymax — 1 and to expand d'lp^d^'^ I dy 
as 

T J (2m) N 

^^ = ^a„(r„+2(x)-r„(x)) (B15) 

rather than improve the numerical conditioning. The expansion (IB 151) explicitly 

builds in the boundary conditions that /dy vanishes at y = 0 and y = ymax and, 

on application of the further condition that ijj^/^\ymax) = 0, means that the expansion 
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X 10“ 



Figure 14. Azimuthal velocities for the (scaled) boundary layer mean flow at 0{mS^/ae^): 
V := /dy — dip^^_j^/dr\r^a^ (thick solid blue line) and the 2nd harmonic - /dy 

(real/imaginary part solid/dashed red line) - plotted against y, the boundary layer variable. 


can be integrated to give 


N 

^ 20 ""^ = anQnjx) where 0„(a;) 

n—Q 


4 _ 3Ti{x) , T 3 (a) 

3 2 6 

3 _ T 2 CX) Ti{x) 

f-T„+i%) T^+3(x)-1 

n+1 ' 2(n+3) 


T„-iix)-l 

2(n-l) 


n = 0 
n = 1 
n > 2 


(B16) 

Using N = 100 gives over 15 decades of drop off in |a„| and the solution is independent 
of ymax well before the actual value ymax = 10 chosen. 


B.3. Solvability condition at 0{S^/e'^) 


The problem for '03™^ 




dy"^ 


= —{Jo + J 2 ) + fl2 


dy^ 


(72- 




10 


dy'^ 


(B17) 
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where the nonlinear terms are 


J2 ■= i- 


■ dipw""'’ (Pip 


(m)* 

10 


dy dy'^ 


_ . # 20 ^ ^ 

dy^ * dy dy'^ 


+ c.c. 




- PPy $° 


dy^ 

(m)* 


y32i2m) 

(0)1^ + - 


(B18) 


with boundary conditions that dip^'/dy(0) = 0 together with V's™ 1 dip^’/dy, and 
d^ip^'^ /dy'^ all vanishing as y —?► 00 . The operator on the LHS of (IB 171) is non-self adjoint 
so we need to find the appropriate adjoint operator to develop a solvability condition. 
Defining the operator 


£1 := 




d 

dy 


(B19) 


then if 4^ := dip^'^/dy (and subscripts indicate derivatives) 


dy = 


$4-^^ - + <S>yy-i> + N{<^y'i> - - N {aQ + 2i2/)$T 


■^c\ ^dy (B 20) 


where 


d^^ nr d'^^ d 


(do -I- 2 * 2 /)$ 


(B21) 


We need now to generate a solution $ of £|$ = 0 such that all the boundary terms 
vanish. The required solution (unique up to arbitrary renormalisation) is, via variation 
of parameters, 


$ := e~2^-yBi[C{y)] 


e5^‘="A*[C(a:)] dx - Ai[C{y)] C e5^‘="Bi[C(x)] dx 


where 


ay) := + 7 ^ + (2*)^/^iVV3,. 


(B22) 


(B 23) 


(2*)2/3 4(2i)2/3 

This has $(0) = $y(0) = 0 and $ ^ 1/// as 2 / —t 00 which ensures all the boundary terms 
in (IB 201) vanish (recall $(0) = dip^^ /dyifi) = 0) and therefore means 


/ $£i^^d2/ = 0. 

Jo dy 

The solvability condition on (IB 171) is then 

/i (72 -I- /2 1)2 = 13 + 14 

where 


T := 




I'l := 


$r_no_d„ 

d2/3 


I 3 ■= 


(B24) 


(B25) 


/ ^Jody, h:= / ^J2dy. 
0 do 
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Since the frequency shift <72 is an imaginary number and the radial flow adjustment is 
real, 


Mih + h]Il) 
Mhll) 


8.47 X 10'‘ > 0 


(B 26) 


as computations give Ii = —24.17—75.26*, I 2 = 140.84 + 20.98*, I 3 = (3.262 + 4.1265*) x 
10® and I 4 = (2.949 — 0.5123*) x 10® (the contribution from the mean flow is an order 
of magnitude larger than that from the 2nd harmonic, both in the same sense). So the 
bifurcation is supercritical with the bifurcating solution branch existing at radial flows 
larger than the critical value. 
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